Parameter identification example

Here is a simple toy model that we use to demonstrate the working of the inference package

$\emptyset \xrightarrow[]{k_1} X \; \; \; \; X \xrightarrow[]{d_1} \emptyset$

Run the MCMC algorithm to identify parameters from the experimental data



In [1]:

    
%matplotlib inline
import bioscrape as bs
from bioscrape.types import Model
from bioscrape.simulator import py_simulate_model
from bioscrape.pid_interfaces import *
from bioscrape.emcee_interface import *
import numpy as np
import pylab as plt
import pandas as pd

pid = initialize_mcmc()
# Import bioscrape XML / bioscrape model object M / SBML file
# M  = bs.types.read_model_from_sbml(filename)
pid.M = Model('models/toy_model.xml')


# Import data from CSV
# Import a CSV file for each experiment run
pid.exp_data = pd.read_csv('test_data.csv', delimiter = '\t', names = ['X','time'], skiprows = 1)
# data_list.append(list(data.get('X')))
# pid.exp_data = np.array(data_list)
pid.measurements = ['X']
pid.time_column = ['time']
pid.prior = {'k1' : ['gaussian', 10, 100, 1e-7],'d1' : ['gaussian', 0.2, 1, 1e-5]}
pid.params_to_estimate = ['k1', 'd1']
pid.nwalkers = 20
pid.init_seed = 0.15
pid.type = 'deterministic'


fitted_model, params = pid.run_mcmc(plot_show = True)









    



C:\ProgramData\Anaconda3\lib\site-packages\emcee\moves\red_blue.py:99: RuntimeWarning: invalid value encountered in double_scalars
  lnpdiff = f + nlp - state.log_prob[j]






    



Successfully completed MCMC parameter identification procedure. Parameter distribution data written to mcmc_results.csv file






    



---------------------------------------------------------------------------
UnboundLocalError                         Traceback (most recent call last)
<ipython-input-1-6af41f3e41eb> in <module>
     29 
     30 
---> 31 fitted_model, params = pid.run_mcmc(plot_show = True)

C:\ProgramData\Anaconda3\lib\site-packages\bioscrape\emcee_interface.py in run_mcmc(self, **kwargs)
     46             self.initial_conditions = self.M.get_species_dictionary()
     47         self.prepare_mcmc(**kwargs)
---> 48         fitted_model, params = self.run_emcee(**kwargs)
     49         return fitted_model.M, params
     50 

C:\ProgramData\Anaconda3\lib\site-packages\bioscrape\emcee_interface.py in run_emcee(self, **kwargs)
    221             f.close()
    222         print('Successfully completed MCMC parameter identification procedure. Parameter distribution data written to mcmc_results.csv file')
--> 223         fitted_model, params = self.plot_mcmc_results(sampler, plot_show)
    224         return fitted_model, params
    225 

C:\ProgramData\Anaconda3\lib\site-packages\bioscrape\emcee_interface.py in plot_mcmc_results(self, sampler, plot_show)
    243             best_p_ind = np.where(n == np.max(n))
    244             if np.shape(best_p_ind)[-1] != 1:
--> 245                 warnings.warn('Multiple parameter values for {0} found with max distribution, choosing the first one. The results might be misleadig.'.format(params_names[i]))
    246                 best_p_ind = np.array(best_p_ind[0].tolist()[0])
    247             assert len(best_p_ind) == 1

UnboundLocalError: local variable 'params_names' referenced before assignment

Check mcmc_results.csv for the results of the MCMC procedure and perform your own analysis.

You can also plot the results as follows



In [ ]:

    
timepoints = pid.timepoints
M_original = Model('models/toy_model.xml')
sim_x = py_simulate_model(timepoints, Model = M_original, stochastic = False)
sim_x_inf = py_simulate_model(timepoints, Model = fitted_model, stochastic = False)
plt.plot(timepoints, sim_x['X'], 'r--', label = 'Original model')
plt.plot(timepoints, pid.exp_data.get('X'), label = 'Experimental data')
plt.plot(timepoints, sim_x_inf['X'], 'k', label = 'Identified model')

plt.legend()

Alll methods above have other advanced options that you can use. Refer to Parameter Identification Tools and Advanced Examples notebook for more details. There are many other tools available such as for multiple initial conditions and timepoints for each trajectory, options for the estimator etc.



In [ ]:

    
fitted_model.get_parameter_values()



In [ ]: